home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dgeev.f < prev    next >
Encoding:
Text File  |  1996-07-19  |  13.6 KB  |  403 lines
  1.       SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
  2.      $                  LDVR, WORK, LWORK, INFO )
  3. *
  4. *  -- LAPACK driver routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       CHARACTER          JOBVL, JOBVR
  11.       INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
  12. *     ..
  13. *     .. Array Arguments ..
  14.       DOUBLE PRECISION   A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
  15.      $                   WI( * ), WORK( * ), WR( * )
  16. *     ..
  17. *
  18. *  Purpose
  19. *  =======
  20. *
  21. *  DGEEV computes for an N-by-N real nonsymmetric matrix A, the
  22. *  eigenvalues and, optionally, the left and/or right eigenvectors.
  23. *
  24. *  The right eigenvector v(j) of A satisfies
  25. *                   A * v(j) = lambda(j) * v(j)
  26. *  where lambda(j) is its eigenvalue.
  27. *  The left eigenvector u(j) of A satisfies
  28. *                u(j)**H * A = lambda(j) * u(j)**H
  29. *  where u(j)**H denotes the conjugate transpose of u(j).
  30. *
  31. *  The computed eigenvectors are normalized to have Euclidean norm
  32. *  equal to 1 and largest component real.
  33. *
  34. *  Arguments
  35. *  =========
  36. *
  37. *  JOBVL   (input) CHARACTER*1
  38. *          = 'N': left eigenvectors of A are not computed;
  39. *          = 'V': left eigenvectors of A are computed.
  40. *
  41. *  JOBVR   (input) CHARACTER*1
  42. *          = 'N': right eigenvectors of A are not computed;
  43. *          = 'V': right eigenvectors of A are computed.
  44. *
  45. *  N       (input) INTEGER
  46. *          The order of the matrix A. N >= 0.
  47. *
  48. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  49. *          On entry, the N-by-N matrix A.
  50. *          On exit, A has been overwritten.
  51. *
  52. *  LDA     (input) INTEGER
  53. *          The leading dimension of the array A.  LDA >= max(1,N).
  54. *
  55. *  WR      (output) DOUBLE PRECISION array, dimension (N)
  56. *  WI      (output) DOUBLE PRECISION array, dimension (N)
  57. *          WR and WI contain the real and imaginary parts,
  58. *          respectively, of the computed eigenvalues.  Complex
  59. *          conjugate pairs of eigenvalues appear consecutively
  60. *          with the eigenvalue having the positive imaginary part
  61. *          first.
  62. *
  63. *  VL      (output) DOUBLE PRECISION array, dimension (LDVL,N)
  64. *          If JOBVL = 'V', the left eigenvectors u(j) are stored one
  65. *          after another in the columns of VL, in the same order
  66. *          as their eigenvalues.
  67. *          If JOBVL = 'N', VL is not referenced.
  68. *          If the j-th eigenvalue is real, then u(j) = VL(:,j),
  69. *          the j-th column of VL.
  70. *          If the j-th and (j+1)-st eigenvalues form a complex
  71. *          conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
  72. *          u(j+1) = VL(:,j) - i*VL(:,j+1).
  73. *
  74. *  LDVL    (input) INTEGER
  75. *          The leading dimension of the array VL.  LDVL >= 1; if
  76. *          JOBVL = 'V', LDVL >= N.
  77. *
  78. *  VR      (output) DOUBLE PRECISION array, dimension (LDVR,N)
  79. *          If JOBVR = 'V', the right eigenvectors v(j) are stored one
  80. *          after another in the columns of VR, in the same order
  81. *          as their eigenvalues.
  82. *          If JOBVR = 'N', VR is not referenced.
  83. *          If the j-th eigenvalue is real, then v(j) = VR(:,j),
  84. *          the j-th column of VR.
  85. *          If the j-th and (j+1)-st eigenvalues form a complex
  86. *          conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
  87. *          v(j+1) = VR(:,j) - i*VR(:,j+1).
  88. *
  89. *  LDVR    (input) INTEGER
  90. *          The leading dimension of the array VR.  LDVR >= 1; if
  91. *          JOBVR = 'V', LDVR >= N.
  92. *
  93. *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
  94. *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  95. *
  96. *  LWORK   (input) INTEGER
  97. *          The dimension of the array WORK.  LWORK >= max(1,3*N), and
  98. *          if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N.  For good
  99. *          performance, LWORK must generally be larger.
  100. *
  101. *  INFO    (output) INTEGER
  102. *          = 0:  successful exit
  103. *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  104. *          > 0:  if INFO = i, the QR algorithm failed to compute all the
  105. *                eigenvalues, and no eigenvectors have been computed;
  106. *                elements i+1:N of WR and WI contain eigenvalues which
  107. *                have converged.
  108. *
  109. *  =====================================================================
  110. *
  111. *     .. Parameters ..
  112.       DOUBLE PRECISION   ZERO, ONE
  113.       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  114. *     ..
  115. *     .. Local Scalars ..
  116.       LOGICAL            SCALEA, WANTVL, WANTVR
  117.       CHARACTER          SIDE
  118.       INTEGER            HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
  119.      $                   MAXB, MAXWRK, MINWRK, NOUT
  120.       DOUBLE PRECISION   ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
  121.      $                   SN
  122. *     ..
  123. *     .. Local Arrays ..
  124.       LOGICAL            SELECT( 1 )
  125.       DOUBLE PRECISION   DUM( 1 )
  126. *     ..
  127. *     .. External Subroutines ..
  128.       EXTERNAL           DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY,
  129.      $                   DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC,
  130.      $                   XERBLA
  131. *     ..
  132. *     .. External Functions ..
  133.       LOGICAL            LSAME
  134.       INTEGER            IDAMAX, ILAENV
  135.       DOUBLE PRECISION   DLAMCH, DLANGE, DLAPY2, DNRM2
  136.       EXTERNAL           LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2,
  137.      $                   DNRM2
  138. *     ..
  139. *     .. Intrinsic Functions ..
  140.       INTRINSIC          MAX, MIN, SQRT
  141. *     ..
  142. *     .. Executable Statements ..
  143. *
  144. *     Test the input arguments
  145. *
  146.       INFO = 0
  147.       WANTVL = LSAME( JOBVL, 'V' )
  148.       WANTVR = LSAME( JOBVR, 'V' )
  149.       IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
  150.          INFO = -1
  151.       ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
  152.          INFO = -2
  153.       ELSE IF( N.LT.0 ) THEN
  154.          INFO = -3
  155.       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  156.          INFO = -5
  157.       ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
  158.          INFO = -9
  159.       ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
  160.          INFO = -11
  161.       END IF
  162. *
  163. *     Compute workspace
  164. *      (Note: Comments in the code beginning "Workspace:" describe the
  165. *       minimal amount of workspace needed at that point in the code,
  166. *       as well as the preferred amount for good performance.
  167. *       NB refers to the optimal block size for the immediately
  168. *       following subroutine, as returned by ILAENV.
  169. *       HSWORK refers to the workspace preferred by DHSEQR, as
  170. *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
  171. *       the worst case.)
  172. *
  173.       MINWRK = 1
  174.       IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
  175.          MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
  176.          IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN
  177.             MINWRK = MAX( 1, 3*N )
  178.             MAXB = MAX( ILAENV( 8, 'DHSEQR', 'EN', N, 1, N, -1 ), 2 )
  179.             K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'DHSEQR', 'EN', N, 1,
  180.      $          N, -1 ) ) )
  181.             HSWORK = MAX( K*( K+2 ), 2*N )
  182.             MAXWRK = MAX( MAXWRK, N+1, N+HSWORK )
  183.          ELSE
  184.             MINWRK = MAX( 1, 4*N )
  185.             MAXWRK = MAX( MAXWRK, 2*N+( N-1 )*
  186.      $               ILAENV( 1, 'DORGHR', ' ', N, 1, N, -1 ) )
  187.             MAXB = MAX( ILAENV( 8, 'DHSEQR', 'SV', N, 1, N, -1 ), 2 )
  188.             K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'DHSEQR', 'SV', N, 1,
  189.      $          N, -1 ) ) )
  190.             HSWORK = MAX( K*( K+2 ), 2*N )
  191.             MAXWRK = MAX( MAXWRK, N+1, N+HSWORK )
  192.             MAXWRK = MAX( MAXWRK, 4*N )
  193.          END IF
  194.          WORK( 1 ) = MAXWRK
  195.       END IF
  196.       IF( LWORK.LT.MINWRK ) THEN
  197.          INFO = -13
  198.       END IF
  199.       IF( INFO.NE.0 ) THEN
  200.          CALL XERBLA( 'DGEEV ', -INFO )
  201.          RETURN
  202.       END IF
  203. *
  204. *     Quick return if possible
  205. *
  206.       IF( N.EQ.0 )
  207.      $   RETURN
  208. *
  209. *     Get machine constants
  210. *
  211.       EPS = DLAMCH( 'P' )
  212.       SMLNUM = DLAMCH( 'S' )
  213.       BIGNUM = ONE / SMLNUM
  214.       CALL DLABAD( SMLNUM, BIGNUM )
  215.       SMLNUM = SQRT( SMLNUM ) / EPS
  216.       BIGNUM = ONE / SMLNUM
  217. *
  218. *     Scale A if max element outside range [SMLNUM,BIGNUM]
  219. *
  220.       ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
  221.       SCALEA = .FALSE.
  222.       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
  223.          SCALEA = .TRUE.
  224.          CSCALE = SMLNUM
  225.       ELSE IF( ANRM.GT.BIGNUM ) THEN
  226.          SCALEA = .TRUE.
  227.          CSCALE = BIGNUM
  228.       END IF
  229.       IF( SCALEA )
  230.      $   CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
  231. *
  232. *     Balance the matrix
  233. *     (Workspace: need N)
  234. *
  235.       IBAL = 1
  236.       CALL DGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
  237. *
  238. *     Reduce to upper Hessenberg form
  239. *     (Workspace: need 3*N, prefer 2*N+N*NB)
  240. *
  241.       ITAU = IBAL + N
  242.       IWRK = ITAU + N
  243.       CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
  244.      $             LWORK-IWRK+1, IERR )
  245. *
  246.       IF( WANTVL ) THEN
  247. *
  248. *        Want left eigenvectors
  249. *        Copy Householder vectors to VL
  250. *
  251.          SIDE = 'L'
  252.          CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL )
  253. *
  254. *        Generate orthogonal matrix in VL
  255. *        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
  256. *
  257.          CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
  258.      $                LWORK-IWRK+1, IERR )
  259. *
  260. *        Perform QR iteration, accumulating Schur vectors in VL
  261. *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
  262. *
  263.          IWRK = ITAU
  264.          CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
  265.      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
  266. *
  267.          IF( WANTVR ) THEN
  268. *
  269. *           Want left and right eigenvectors
  270. *           Copy Schur vectors to VR
  271. *
  272.             SIDE = 'B'
  273.             CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
  274.          END IF
  275. *
  276.       ELSE IF( WANTVR ) THEN
  277. *
  278. *        Want right eigenvectors
  279. *        Copy Householder vectors to VR
  280. *
  281.          SIDE = 'R'
  282.          CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR )
  283. *
  284. *        Generate orthogonal matrix in VR
  285. *        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
  286. *
  287.          CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
  288.      $                LWORK-IWRK+1, IERR )
  289. *
  290. *        Perform QR iteration, accumulating Schur vectors in VR
  291. *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
  292. *
  293.          IWRK = ITAU
  294.          CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
  295.      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
  296. *
  297.       ELSE
  298. *
  299. *        Compute eigenvalues only
  300. *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
  301. *
  302.          IWRK = ITAU
  303.          CALL DHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
  304.      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
  305.       END IF
  306. *
  307. *     If INFO > 0 from DHSEQR, then quit
  308. *
  309.       IF( INFO.GT.0 )
  310.      $   GO TO 50
  311. *
  312.       IF( WANTVL .OR. WANTVR ) THEN
  313. *
  314. *        Compute left and/or right eigenvectors
  315. *        (Workspace: need 4*N)
  316. *
  317.          CALL DTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
  318.      $                N, NOUT, WORK( IWRK ), IERR )
  319.       END IF
  320. *
  321.       IF( WANTVL ) THEN
  322. *
  323. *        Undo balancing of left eigenvectors
  324. *        (Workspace: need N)
  325. *
  326.          CALL DGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL,
  327.      $                IERR )
  328. *
  329. *        Normalize left eigenvectors and make largest component real
  330. *
  331.          DO 20 I = 1, N
  332.             IF( WI( I ).EQ.ZERO ) THEN
  333.                SCL = ONE / DNRM2( N, VL( 1, I ), 1 )
  334.                CALL DSCAL( N, SCL, VL( 1, I ), 1 )
  335.             ELSE IF( WI( I ).GT.ZERO ) THEN
  336.                SCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ),
  337.      $               DNRM2( N, VL( 1, I+1 ), 1 ) )
  338.                CALL DSCAL( N, SCL, VL( 1, I ), 1 )
  339.                CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 )
  340.                DO 10 K = 1, N
  341.                   WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2
  342.    10          CONTINUE
  343.                K = IDAMAX( N, WORK( IWRK ), 1 )
  344.                CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R )
  345.                CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
  346.                VL( K, I+1 ) = ZERO
  347.             END IF
  348.    20    CONTINUE
  349.       END IF
  350. *
  351.       IF( WANTVR ) THEN
  352. *
  353. *        Undo balancing of right eigenvectors
  354. *        (Workspace: need N)
  355. *
  356.          CALL DGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR,
  357.      $                IERR )
  358. *
  359. *        Normalize right eigenvectors and make largest component real
  360. *
  361.          DO 40 I = 1, N
  362.             IF( WI( I ).EQ.ZERO ) THEN
  363.                SCL = ONE / DNRM2( N, VR( 1, I ), 1 )
  364.                CALL DSCAL( N, SCL, VR( 1, I ), 1 )
  365.             ELSE IF( WI( I ).GT.ZERO ) THEN
  366.                SCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ),
  367.      $               DNRM2( N, VR( 1, I+1 ), 1 ) )
  368.                CALL DSCAL( N, SCL, VR( 1, I ), 1 )
  369.                CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 )
  370.                DO 30 K = 1, N
  371.                   WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2
  372.    30          CONTINUE
  373.                K = IDAMAX( N, WORK( IWRK ), 1 )
  374.                CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R )
  375.                CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
  376.                VR( K, I+1 ) = ZERO
  377.             END IF
  378.    40    CONTINUE
  379.       END IF
  380. *
  381. *     Undo scaling if necessary
  382. *
  383.    50 CONTINUE
  384.       IF( SCALEA ) THEN
  385.          CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
  386.      $                MAX( N-INFO, 1 ), IERR )
  387.          CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
  388.      $                MAX( N-INFO, 1 ), IERR )
  389.          IF( INFO.GT.0 ) THEN
  390.             CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
  391.      $                   IERR )
  392.             CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
  393.      $                   IERR )
  394.          END IF
  395.       END IF
  396. *
  397.       WORK( 1 ) = MAXWRK
  398.       RETURN
  399. *
  400. *     End of DGEEV
  401. *
  402.       END
  403.